MASTER  COPY 


KEEP  FOR  REPRODUCTION  PURPOSES 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  NO.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comment  regarding  this  burden  estimates  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188),  Washington,  DC  20503. 


1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

February  28,  1997  Final,  1994  -  1996 


4.  TITLE  AND  SUBTITLE  5.  FUNDING  NUMBERS 

A  Study  of  Planetary  Boundary  Layer  Forcing  on  the  DAAH04-94-G-0049 

Turbulent  Flow  Characteristics  of  Tree  Canopies  Using 


6.  AUTHOR(S) 

David  R.  Miller 
X.  Harrison  Yang 


7.  PERFORMING  ORGANIZATION  NAMES(S)  AND  ADDRESS(ES) 

University  of  Connecticut 

Natural  Resources  Management  &  Engineering  Department 
1376  Storrs  Road 
Storrs ,  CT  06269-4087 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

U.S.  Army  Research  Office 
P.O.  Box  12211 

Research  Triangle  Park,,  NC  27709-2211 


10.  SPONSORING /MONITORING 
AGENCY  REPORT  NUMBER 


11.  SUPPLEMENTARY  NOTES 


The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  authors)  and  should  not  be  construed  as 
an  official  Department  of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12  b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution  unlimited. 


13.  ABSTRACT  (Maximum  200  words) 


19970321  028 


A  non-local  closure  numerical  model  based  on  the  Transilient  Turbulence  Theory  (TTT)  (Stull,  1984-1986)  has  been 
developed  and  verified  in  and  above  a  homogeneous  tree  canopy.  The  first  order,  volumetricaily  averaged  governing  equations  for 
momentum,  heat,  and  humidity  were  coupled  dynamically  with  canopy  processes  such  as  evapotranspiration  and  radiation  within 
the  canopy.  The  simulated  profiles  of  momentum,  heat,  and  humidity  were  compared  to  measurements  in  a  homogeneous  orchard 
during  project  WTND.  The  model  provides  evidence  that  non-local  turbulence  closure  by  TTT  is  a  consistent  and  efficient  technique 
to  parameterize  the  turbulent  transport  near  the  ground  in  canopies  as  well  as  above,  and  can  be  incorporated  into  numerical  models. 

The  spatial  resolution  in  surface  parameterizations  for  spatially  distributed  models  was  quantified  using  spatial 
autocorrelations.  The  resulting  error  prediction  model  was  used  to  quantify  the  statistical  errors  due  averaging  different  patch  sizes 
and  contrasts.  Sensitivity  of  process  models  to  these  averaging  errors  in  boundary  conditions  was  quantified  using  probable  error 
vectors. 

Overall  these  efforts  have  rigorously  set  the  stage  for  future  use  of  the  non-local  turbulent  transport  approach  and  as  an 
efficient  and  accurate  method  to  model  spatially  distributed  canopy  flow  in  three  dimensions 


DTIC  QUALITY  INSPECTED  2 


14.  SUBJECT  TERMS 

Transilient  Turbulence,  Non-local  Closure,  Tree  Canopy,  Flow 
Modelling,  Spatial  Averaging 


15.  NUMBER  IF  PAGES 

31 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION 
OR  REPORT 

UNCLASSIFIED 


18.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

UNCLASSIFIED 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

UNCLASSIFIED 


20.  LIMITATION  OF  ABSTRACT 


0 


A  Study  of  Planetary  Boundary  Layer  Forcing  on  Turbulent  Flow 
Characteristics  in  Tree  Canopies  Using  Project  Wind  Data 

Final  Progress  Report 


David  R.  Miller  and  X.  Harrison  Yang 


February  28, 1997 


U.  S.  Army  Research  Office 
Grant  Number  DAAH04-94-G-0049 


Department  of  Natural  Resources  Management  and  Engineering 
University  of  Connecticut 


APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED. 


THE  VIEWS,  OPINIONS,  AND  /OR  FINDINGS  CONTAINED  IN  THIS  REPORT 
ARE  THOSE  OF  THE  AUTHORS(S)  AND  SHOULD  NOT  BE  CONSTRUED  AS  AN 
OFFICIAL  DEPARTMENT  OF  THE  ARMY  POSITION,  POLICY,  OR  DECISION, 
UNLESS  SO  DESIGNATED  BY  OTHER  DOCUMENTATION. 


1 


r 


A.  STATEMENT  OF  THE  PROBLEM  STUDIED. 

The  primary  purpose  of  this  project  was  to  develop  a  Transilient  Turbulence  non¬ 
local  closure  model  of  the  wind  and  energy  flux  in  tree  canopies  which  would  realistically 
couple  the  roughness  layer  wind  field  with  the  wind  field  above  the  canopy. 


B.  SUMMARY  OF  MOST  IMPORTANT  RESULTS 

A  non-local  closure  numerical  model  based  on  the  Transilient  Turbulence  Theory 
(TTT)  (Stull,  1984-1986)  has  been  developed  and  verified  in  and  above  a  homogeneous 
tree  canopy.  The  first-order,  volumetrically  averaged  governing  equations  for 
momentum,  heat,  and  humidity  were  coupled  dynamically  with  canopy  processes  such  as 
evapotranspiration  and  radiation  within  the  canopy.  The  simulated  profiles  of 
momentum,  heat,  and  humidity  were  compared  to  measurements  in  a  homogeneous 
orchard  during  project  WIND.  The  model  provides  evidence  that  nonlocal  turbulence 
closure  by  TTT  is  a  consistent  and  efficient  technique  to  parameterize  the  turbulent 
transport  near  the  ground  in  canopies  as  well  as  above,  and  can  be  incorporated  into 
numerical  models.  The  Appendix  of  this  report  contains  the  documentation  of  the 
model. 

The  project  WIND  data  sets  from  the  almond  orchard  were  summarized  into  wind 
direction  and  stability  classes  and  utilized  to  verify  the  model. 

The  TTT  approach  to  modeling  canopy  fluxes  and  coupling  them  to  the  PBL 
turned  out  to  be  very  sensitive  to  accurate  specification  of  the  surface  boundary 
conditions.  Therefore  an  investigation  into  the  accuracy  of  spatially  averaging  the  surface 
boundary  conditions  was  conducted.  The  spatial  resolution  in  surface  parameterizations 
for  spatially  distributed  models  was  quantified  using  spatial  autocorrelations.  The 
resulting  error  prediction  model  was  used  to  quantify  the  statistical  errors  due  averaging 
different  patch  sizes  and  contrasts.  Sensitivity  of  process  models  to  these  averaging 
errors  in  boundary  conditions  was  quantified  using  probable  error  vectors. 

Overall  these  efforts  have  rigorously  set  the  stage  for  future  use  of  the  non-  local 
turbulent  transport  approach  and  as  an  efficient  and  accurate  method  to  model  spatially 
distributed  canopy  flow  in  three  dimensions. 
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A.  Physical  Configuration  of  the  Model 


The  model  is  one  dimensional  and  horizontally  homogeneous  with  different  vertical  spacings  and  timesteps 
in  three  regions.  The  first  is  called  the  roughness  sublayer  region  which  extends  from  the  ground  surface  up 
to  four  tree  heights.  In  this  region,  the  vertical  spacing  is  Az  =  1  m  and  the  timestep  is  At  =  3.6  sec . 

The  second  region  extend  from  four  to  ten  tree  heights.  The  vertical  spacing  is  Az  =  10  m ,  and  the 
timestep  is  At  =  36  sec .  The  combination  of  the  first  and  second  regions  may  be  seen  as  the  surface 
layer. 

From  the  ten  tree  heights  up  to  the  top  of  the  planetary  boundary  layer  (PBL:  ~3  km),  the  vertical 

spacing  is  Az  =  1 00  m ,  and  the  timestep  is  At  =  360  sec  . 

The  simulation  area  is  divided  into  three  regions  because  the  distinguishing  characteristics  of  these 
layers  are  different  in  response  to  the  thermodynamic  forcings  within  the  PBL.  Within  the  roughness 
sublayer,  the  turbulence  is  often  characterized  by  high  intensity  and  large  skewness  and  kurtosis,  which 
indicates  that  the  turbulence  is  non-Gaussian  and  intermittent.  This  is  true  in  all  sizes  of  canopies  from  tall 
forests  (Baldocchi  and  Hutchison,  1987;  Baldocchi  and  Meyers,  1988)  to  short  grass  (Aylor  et  al. ,  1993). 

Above  the  roughness  sublayer,  however,  the  turbulence  can  be  properly  described  by  the  Gaussian 
distribution  as  over  a  flat  plain  (Kaimal  and  Finnigan,  1994;  Seginer  et  al. ,  1976). 

The  use  of  an  intermediate  region  is  due  to  both  physical  and  numerical  concerns.  The  surface  layer, 
which  extends  from  the  ground  surface  to  about  100  m,  is  observed  to  be  decoupled  from  the  mixing  layer 
at  night  and  only  strong  instability  can  break  the  surface  inversion  layer  produced  at  night,  requiring 
special  consideration  in  the  model.  Also  in  order  to  keep  the  numerical  scheme  smooth  in  terms  of  vertical 
spacing  and  timestep,  it  is  necessary  to  have  this  intermediate  region  to  avoid  jumping  between  two 
extreme  timesteps  (3.6  sec  and  360  sec). 

Reference 

Aylor,  D.  E.,  Wang,  Y.  Y.,  Miller,  D.  R.:  1993,  ‘Intermittent  wind  close  to  the  ground  within  a  grass 
canopy’,  Bound-Layer  Meteorol.  66,  427-448. 
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Baldocchi,  D.  D.,  and  Hutchison,  B.  A.:  1987,  ‘Turbulence  in  an  almond  orchard:  Vertical  variation  in 
turbulent  statistics’,  Bound-Layer  Meteor ol.  40,  127-146. 

Baldocchi,  D.  D.,  and  Meyers,  T.  P.:  1988,  ‘Turbulence  structure  in  a  deciduous  forest5,  Bound.-Layer 
Meteorol.  43,  345-364. 

Kaimal,  J.  C.  and  Finnigan,  J.  L.:  1994,  Atmospheric  Boundary  Layer  Flows  -  Their  Structure  and 
Measurement ,  Oxford  University  Press,  New  York,  289  pp. 

Seginer,  I.,  Mulheam,  P.  J.,  Bradley,  E.  F.,  and  Finnigan,  J.  J.:  1976,  Turbulent  flow  in  a  model  plant 
canopy5,  Bound.-Layer  Meteorol.  10,  423-453. 
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Flowchart  of  the  Model' 


*  This  is  the  flowchart  for  the  main  (control)  program  of  the  model,  which  is  named  as  CAN_MAIN.  The 
programs  and  symbols  are  explained  in  the  following  text. 


9 


Leaf  Area  Index  and  Leaf  Area  Density  (LF) 


Algorithm 

The  cumulative  leaf  area  index  ( L :  m2  m'2)  of  a  canopy  is  calculated  using  Weibull’s  cumulative 
distribution  function  following  Yang  et  al  (1995), 


L(z )  =  LAIO  x  <  1  -  exp 


(LF.l) 


where  LAI  represents  the  leaf  area  index  excluding  stems,  z  is  the  height  in  meters  above  the  ground,  H  is 
the  tree  height  in  meters,  and  b  and  c  are  the  two  parameters  specific  to  the  canopy  structures.  The  program 
for  Eq.(LF.l)  is  named  as  CUMLAI(z). 

Eq.(LF.l)  may  be  properly  used  for  the  distribution  of  radiation.  In  consideration  of  the  drag  of  the 
canopy,  stems  and  trunks  should  be  also  included,  although  they  do  not  participate  in  the  radiation  budget. 
By  adding  one  linear  term  for  the  trunk,  Eq.(LF.l)  may  be  modified  as 


L(z)  =  LAI1  x 


+  0..3-(l  -z/H) 


(LF.2) 


where  LAI1  is  leaf  area  index  including  stems.  The  program  is  named  as  CUMLAIO(z). 
The  leaf  area  density  ( a :  m2  m'3)  is  evaluated  according  to  its  definition 


a(z)  = 


dL(z ) 
dz 


which  yields,  from  the  discretization  of  equation  (LF.2), 


a(z) 


L(z  +  Az/2)  -  L(z  -  Az/2) 
A z 


(LF.3) 


(LF.4) 


where  Az  is  the  vertical  spacing.  The  corresponding  programs  of  leaf  area  density  to  Eqs.(LF.l)  and  (LF.2) 
are  named 


and 


LEAFAREADENSITY(z) 
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LEAFAREADENSITYO(z), 


respectively. 


Reference 

Yang,  X.,  Miller,  D.  R.,  and  Witcosky,  J.  J.:  1994,  ‘Characterization  of  hardwood  forest  canopies  in  the 
eastern  United  States’,  21s *  Conference  on  Agriculture  and  Forest  Meteorology,  March  7-11, 
1994,  130-133. 
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Boundary  Layer  Resistance  (BR) 


Algorithm 


Following  Goudriaan  (1977),  the  boundary  layer  resistance  (rfr.  s  m'1)  of  the  foliage  at  each  grid  box  may 
be  estimated  as 


,  dl 

rb  =  0.5  x  1.8  x  102  x 


(BR.1) 


where  dl  is  the  characteristic  width  of  leaves,  and  M  = 


-\/u2  +  v2 


is  the  mean  wind  speed. 


Program 


The  function  boundary  layer  resistance  is  named  as 
BL_RESISTANCE(u,v) 

where  u  and  v  represent  the  two  components  of  the  wind  velocity. 


Reference 

Goudriaan,  1977:  Crop  Micrometeorology:  A  Simulation  Study,  249  pp. 
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Stomata!  Resistance  (SR) 


Algorithm 


As  has  been  discussed  by  Norman  (1979),  the  stomatal  resistance  (rs)  strongly  depends  on  the  irradiance 
on  the  leaves  and  the  bio-physical  conditions  of  the  leaves.  It  may  be  expressed  as 


r^MKVMMUv,.)  (SR.D 

where  f,f2,  an df3  are  functional  representations,  Rv  is  the  visible  (or  photosynthetically  active)  radiation 
the  leaves  received,  y/lw  ls  the  leaf  water  potential.  Norman  (1979)  suggests  that  the  three  functions  be 
given  as 


- 

(  *  ^ 

- 

K 

-l 

+  Rv 

\  ^s,min  ) 

- 

(  *  > 

- 

r 

A-,min 

K 

-l 

+  Rv 

+  K 

r 

cuticle 

\  ^v,min  ) 

/2(o= 


1  +  m-\ 


1  +  m- 


Op  ~  %,min 
| _^c,U  ®c, min  _ 

%,min  ~  0c 
^c.min  ~  @c,L  _ 


for  min 


for  0C  <  6cmm 


Mvlw)= 


1  Vh,  >  V™ 

V'-  -  j "JL  wm>w  >vm 

(2)  Tlw  ^  Tlw  ^  T  lw 


-  n 


. 0  Vtw  <  Vi 

where  the  parameters  in  the  three  functions  are  listed  in  Table  I. 


(2) 


(SR.2) 


(SR.3) 


(SR.4) 


Table  I 

Parameter  Values  and  Interpretations  in  the  Three  Functions  for  the  Calculation  of  the  Stomatal  Resistance. 
Symbols  Values  Interpretation 
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f(Rv)  (equation  (SR.2)) 


fs,min 

200  s  m'1 

minimum  stomatal  resistance 

K 

100  W  nr2 

critical  visible  radiation 

* 

400  s  m'1 

critical  stomatal  resistance 

rs 

^cuticle 

3000  s  nr1 

cuticular  resistance 

f2{0c)  (equation  (SR.3)) 

m 

5 

fitting  parameter 

^C,min 

30  °C 

temperature  corresponding  to  the  minimal 

MM 

dc,L 

15  °C 

lower  reference  temperature 

9c,U 

45  °C 

upper  reference  temperature 

Mwi w)  (equation  (SR.4)) 

-12  bars 

upper  reference  water  potential 

-17  bars 

lower  reference  water  potential 

Program 


The  subroutine  is  named  as 

STOMA_RESISTANCE(rs,RV,Tc) 

where  rs  represents  the  stomatal  resistance  as  the  output.  The  inputs  are  the  visible  radiation  (RV)  and  the 
canopy  temperature  (Tc). 

Reference 

Norman,  J.  M.:  1979,  ‘Modeling  the  complete  crop  canopy’,  In:  Modification  of  the  Aerial  Environment  of 
Crops ,  ASAE  Monograph ,  249-3 1 1 . 
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Irradiance  (IR) 


Algorithms 


Following  Goudriaan  (1977),  the  net  irradiance  within  each  grid  box  is  represented  as 

^net  =  ^v^v,net  +  ^NIR  ^NIR,net  +  ^L,net  (1^*0 

where  Rnet,  /?vnet,  and  tfNIIMet  represent  the  (total)  net,  the  visible,  and  the  near  infrared  irradiances, 
respectively;  and  Rh  represents  the  net  longwave  radiation  in  the  grid  box;  and  aNIR  represent  the 
absorption  coefficients  of  the  leaves  for  the  visible  and  the  near  infrared  wavebands,  respectively.  The 
ultraviolet  irradiance  is  not  included  in  the  study  for  its  negligible  contribution  to  the  total  radiation 
(Goudriaan,  1977). 

In  equation  (IR.l),  the  visible  (/?vnet)  and  near  infrared  (^NIRinet)  irradiances  are  further  categorized  into 
direct  (or  beam)  and  diffuse  parts,  respectively.  The  profiles  of  direct  and  diffuse  irradiances  within  the 
canopy  are  calculated  following  exponential  extinction  relationship 

R  =  Rsop  expf-  KL)  (IR.2) 

where  R  represents  the  direct  (beam)  or  diffuse  irradiance  of  either  the  visible  or  the  near  infrared,  Rtop 
represents  the  incoming  radiation  at  the  top  of  the  canopy,  K  represents  the  extinction  coefficient  of  the 
radiation  within  the  canopy,  and  L  is  the  cumulative  leaf  area  index  defined  in  equation  (LF.l). 

The  study  is  designed  for  a  clear  sky  so  that  the  visible  and  the  near  infrared  radiation  contribute 
evenly  to  the  total  shortwave  radiation,  given  that  the  ultraviolet  part  is  neglected. 

The  partition  of  the  direct  and  the  diffuse  radiation  at  the  top  of  the  canopy  follows  the  relationship 


0.0837 

/d=——  +  03 


(IR.3) 


fB  =!-//>  <IR-4) 

where  7b  and  fB  are  the  partition  coefficients  for  the  diffuse  and  direct  radiation,  respectively.  Equations 
(IR.3)  and  (IR.4)  were  obtained  by  linear  regression  according  to  the  numerical  values  given  by  Goudriaan 


(1977).  The  solar  altitude,  sina,  is  given  by  Bras  (1990)  as 
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sina  =  sin  Jsin®  +  cost)  cos  <I>  cos 


(IR.5) 


where  S,  0.  and  y/  are  declination  of  the  sun,  local  latitude,  and  the  hour  angle  of  the  sun,  respectively.  The 
declination  of  the  sun  is  given  as 


5  = 


23.45;r 

180 


cos 


2  n 
365 


(172  -D) 


(1R.6) 


where  D  is  the  Julian  day  (1  <  D  <  365  or  366 )  and  8 is  in  radians. 
The  local  hour  angle,  0  <\f/<  2 71 ,  is  given  as 


y/  =  ) 


n 

Y2 

n 

12 


(rs+12-AT) 

(rs-12-A7-) 


when  the  sun  is  east  of  the  observer' s  longitude, 
when  the  sun  is  west  of  the  observer' s  longitude, 


(IR.7) 


where  T$  is  the  standard  time  in  the  time  zone  of  the  observer  in  hours  counted  from  midnight  ( 0.00 - 
23.59),  AT  is  the  time  difference  between  standard  and  local  longitude  in  hours  given  as 


at  =  ^(0s-Ol),  (IR.8) 

where  /  =  -1  for  west  longitude  and  /  =  1  for  east  longitude,  relative  to  Greenwich,  Os  the  standard 
meridian  (meridian  where  the  observer’s  time  zone  is  centered),  and  61  is  the  longitude  of  the  observer 
meridian. 

In  equation  (IR.2),  the  extinction  coefficient  for  the  direct  (beam)  radiation  is  given  as 

Hh£ 

sin  a 

where  O(a)  is  the  average  projection  given  as 

<9(/?)  =  0,  +0.877(l  -20,)  sin  a,  (IR.10) 

with 


\'/2 


(IR.9) 


Ox  =  0.5  -  0.633^  -  033x1  '  0-4  <  XL  <  0-6 , 


(IR.ll) 
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1 

^=±2 


(IR.12) 


3  6  9 

0.1 34  -  X  F{X)  +  0.366  -  X  F(X)  +  0.5  -  X  F(X) 

X=\  X=4  X=7 

where  0]  is  an  intermediate  function,  and  XL  represents  the  deviation  of  the  actual  leaf  angle  distribution 
from  the  spherical  one.  The  sign  of  XL  is  the  same  as  that  of  the  last  tenn. 

The  function  F(X)  in  equation  (IR.12)  describes  the  leaf  inclination  distribution.  For  this  study,  it  is 
assumed  that  the  leaf  inclination  is  isotropic  so  that  F{X)  =  1/9  in  each  leaf  angle  class.  However,  it  may 
be  more  accurate  for  the  calculation  of  XL  if  the  measurement  of  F(X)  is  known. 

In  equation  (IR.9),  cris  defined  from  the  simplification  of  the  leaf  reflection  (p )  and  transmission  (r) 
coefficients  as  p  —  T  —  0.5<T . 

The  extinction  coefficient  for  the  diffuse  radiation  within  the  canopy  is  given  as 

Kd  =  S,(A)exp{-  Kb(p,)L]\ ,  (IIU3) 

where  Bu(Pn)  denotes  the  distribution  of  the  diffuse  light  over  the  inclination  class  J3n  under  the  assumption 
that  the  sky  has  a  uniform  radiance, 

=  -|jcos  Y  ~cos  -  |  for n  =  1,2, •••,9 .  (IR.14) 

The  inclination  of  the  diffuse  light  is  given  as 

^  =  4{l0(”-1)+5}  forn  =  1’2’-’9-  (IR15) 

The  extinction  coefficient  K\y  in  equation  (IR.13)  is  given  by  equation  (IR.9)  except  that  the  solar 
inclination  angle  a  is  now  replaced  by  the  angle  of  the  diffuse  light  of  the  nine  inclination  classes  fin  (n  = 
1,2,...,  9). 

For  the  calculation  of  longwave  radiation,  it  is  assumed  that  the  canopy  as  a  whole  exchanges  energy 
first  with  the  clear  sky  and  then  redistribute  the  gain/loss  of  the  energy  within  the  canopy  according  to 
exponential  extinction  relationship.  Mathematically,  the  exchange  of  longwave  radiation  between  the  sky 


and  the  canopy  is  written  as 


(IR.16) 


BL  =  *sky°3*,  -  £<.°dt  > 

where  ssky  and  ec  are  the  emissivities  of  the  sky  and  the  canopy,  respectively,  cr  is  the  Stefan-Boltzmann 
constant  which  has  the  value  5.67  x  1 0  8  W  m  '2 , 3ky  is  the  apparent  potential  temperature  of  the  sky, 
and  0m  is  the  average  temperature  of  the  canopy. 

According  to  Idso  (1981), 

=0.70  +  5.95x1  O'7  e-exp(l500/6>sky)  (IR.17) 

where  0sky  is  the  potential  air  temperature  in  Kelvin  at  2  m  height  above  the  canopy.  (From  the  test  run  of 
the  study,  equation  (IR.17)  may  be  subject  to  adjustment  in  accordance  with  seasons  and  locations). 

The  redistribution  of  the  longwave  exchange  within  the  canopy  may  be  estimated  as 

R-L  =  Bl  exp(-  KlL),  (IR.18) 

where  K £  is  the  extinction  coefficient  which  has  the  value  0.81  for  the  longwave  radiation  within  the 
canopy  (Goudriaan,  1977). 

Overall,  the  net  irradiance  within  /th  grid  box  in  the  canopy  is  calculated  by  equation  (IR.l).  The  net 
irradiance  Rvneti  tfNIR  net,  and  ^L.net  within  each  grid  box  are  calculated  as 

-  R.w  =  IA5-7)  •  (,R-20) 

where  x  represents  different  wave  bands  (v:  visible;  NIR:  near  infared;  L:  longwave).  The  irradiance  Rxj+  m 
has  been  prescribed  by  equation  (IR.2)  for  the  visible  and  near  infrared  wavebands  of  direct  and  diffuse 
radiation  and  equation  (IR.18)  for  the  longwave  radiation.  The  half  grid  number  indicates  that  the  radiation 
is  computed  at  even  number  of  heights  (0,  2,  4,  6,  8  m)  which  are  presumably  the  boundaries  of  the  grid 
boxes. 

The  total  radiation  reaching  at  the  top  of  the  canopy  may  be  evaluated  by  theoretical  formula  or  by  the 
observation  as  the  input.  It  is  summarized  by  Bras  (1990)  that  the  clear  sky  shortwave  radiation,  Ic,  after 
accounting  for  atmospheric  effects,  as 

Ic  =  — f-(sina)exp(-«a1w),  (IR.21) 
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where  W0  is  the  solar  constant  (=  1353  W  m'2),  r  is  the  ratio  of  actual  earth-sun  distance  to  mean  earth-sun 
distance,  given  by 


r  =  1.0  +  0.017 cos|  |^(1 86  -  D) 
365 


(IR.22) 


n  is  a  turbidity  factor  (=  2  for  clear  mountain  air),  ax  is  the  molecular  scattering  coefficient  defined  as 

ax  =  0.128  -0.0541og10m,  (IR.23) 

and  m  is  the  effective  thickness  of  the  atmosphere  (or  called  optical  air  mass)  estimated  as 


m  =  [sin  a  +  0.1 5(a  +  3.885)“'  253  ]“ 


(IR.24) 


Program 

The  subroutine  to  calculate  the  total  net  irradiance  (equation  (IR.l))  is 
NET_RAD(RN,RV,Ta,Tm,e,j) 

where  RN  is  the  total  net  irradiance,  RV  is  the  net  visible  irradiance  which  is  for  the  calculation  of  the 
stomatal  resistance  (see  equation  (SR.2))  and  subroutine  STOMA_RESISTANCE),  Ta  is  the  potential  air 
temperature  at  2  m  above  the  canopy,  Tm  is  the  average  canopy  temperature,  and  e  is  the  vapor  pressure, 
and  j  is  the  index  of  the  grid  box. 

The  subroutine  to  calculate  the  diffuse  irradiance  (equation  (IR.2))  is 

DIFFUSE(RD,RD_IN, rhotau, j) 

where  RD  is  the  net  output  diffuse  (visible  or  near  infrared)  radiation,  RD  IN  is  the  input  diffuse  (visible 
or  near  infrared)  radiation  at  the  top  of  the  canopy,  rhotau  is  the  summation  of  the  reflection  and 
transmission  coefficients,  and  j  is  the  grid  box  index. 

The  subroutine  to  calculate  the  net  direct  irradiance  (equation  (IR.2))  is 

DIRECT(RB,RB_lN,sina, rhotau, j) 

where  RB  is  the  net  output  direct  (visible  or  near  infrared)  radiation,  RB  IN  is  the  input  diffuse  (visible  or 
near  infrared)  radiation  at  the  top  of  the  canopy,  and  sina  is  the  solar  altitude. 

The  longwave  radiation  (equation  (IR.  1 8))  is  calculated  by  the  subroutine 
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LONGWAVE(RL,Ta,Tm,e,j) 


where  RL  is  the  net  output  longwave  radiation.  The  rest  variables  were  explained  above. 

The  extinction  coefficient  for  the  diffuse  radiation  (equation  (IR.  13))  is  given  in  the  subroutine 

EXTINCTD(Kd,rhotau) 

where  Kd  is  the  output  extinction  coefficient  for  the  diffuse  radiation. 

The  extinction  coefficient  for  the  direct  radiation  (equation  (IR.9))  is  calculated  by  the  subroutine 

EXTINCTB(Kf,sina,rhotau) 

where  Kf  is  the  output  extinction  coefficient  for  the  direct  radiation. 

The  distribution  of  the  diffuse  light  over  the  nine  inclination  class  (equation  (IR.14))  is  given  in  the 
subroutine 

DISJTBL(Bu,i) 

where  Bu  is  the  distribution  function. 

The  solar  radiation  at  the  top  of  the  canopy  (equation  (IR.21))  and  the  partition  of  the  direct  and 
diffuse  of  both  visible  and  near  infrared  radiation  are  given  in  the  subroutine 

SOLAR_RAD(RD_V,RD_NIR,RB_V,RB_NIR,sina) 


Reference 

Bras,  R.  L.:  1990,  Hydrology.  Andison- Wesley  Publishing  Company,  New  York,  643  pp. 

Goudriaan,  J.:  1977,  Crop  Micrometeorology:  A  Simulation  Study.  PUDOC,  Center  for  Agric.  Publ.  and 
Doc.  Wageningen,  The  Netherland,  249  pp. 

Idso,  S.  B.:  1981,  ‘A  set  of  equations  for  full  spectrum  and  8-  to  14-mm  and  10.5-  to  12.5-mm,  thermal 
radiation  from  cloudiness  skies’,  Water  Resources  Res.  17,  295-304. 
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Canopy  Temperature  (CT) 
Algorithm 


Non-steady  State 

According  to  Jones  (1992),  the  leaf  temperature  within  each  grid  box  at  non-steady  state  is  written  as 

dOc  R  -  SH  -  LH 

— ^  =  -  (CT.l) 

dt  PcCJc 

where  Tc  is  the  canopy  temperature,  t  is  the  time,  SH  is  the  sensible  heat,  LH  is  the  latent  heat,  pc  and  Cp  C 
are  the  density  and  specific  heat  capacity,  respectively,  of  leaf  tissue.  lc  is  the  thickness  of  a  flat  leaf  id! 4 
for  a  cylinder  or  d/6  for  a  sphere,  where  d  is  the  diameter),  and  Az  is  the  spatial  interval  in  the  vertical 
direction. 

The  sensible  heat,  SH,  is  given  by  Norman  (1979)  and  Jones  (1992)  as 


SH  = 


2aAzpaCpa 


(CT.2) 


where  pa  and  Cp >a  are  the  density  and  specific  heat,  respectively,  of  the  air.  The  factor  of  two  corresponds 
to  the  two  sides  of  the  leaves. 

The  latent  heat,  LH,  is  given  as 


LH  - 


(CT.3) 


where  s  is  the  ratio  of  the  molecular  weights  of  water  vapor  and  dry  air  (  =  0.622  ),  X  is  the  latent  heat  of 


evaporation  (=2.45x  1 06  J  kg’1),  P  is  the  surface  atmospheric  pressure  (=  105  Pa),  es(6c)  is  the 

saturation  vapor  pressure  at  the  canopy  temperature  Qc,  and  e  is  the  in  situ  vapor  pressure. 

Given  the  specific  humidity  of  the  air  q  and  the  surface  atmospheric  pressure  P,  the  vapor  pressure  in 
equation  (CT.3)  is  calculated  as 


e  = 


P 

0.622  q 


(CT.5) 
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The  saturation  vapor  pressure,  es{0a ),  is  given  as 


*,(<0  =  611 


x  exp 


17.67(3,  -273.15) 

a.  -29.66 


(CT.6) 


Program 

The  main  subroutine  to  calculate  the  canopy  temperature  (equation3  (CT.1-3))  is  named  as 

CANOPY_TMP(thetac,  SH,  LH,  thetav,  RN,  u,  v,  q,  tmp  a,  flxRs,  abs_Rs, 

TMP  RL,  TMP  RV,  M,  Mh) 

where  SH  is  the  sensible  heat,  LH  represents  the  latent  heat,  thetav  represents  the  air  temperature,  RN  is 
the  net  irradiance,  u  and  v  are  the  velocity  components  of  winds,  q  represents  the  specific  humidity,  tmp_a 
is  the  output  of  the  leaf  area  density,  flx_Rs  is  the  output  of  the  shortwave  radiation  flux,  abs_Rs  is  the 
output  of  the  absorption  of  the  shortwave  radiation,  TMP_RL  is  the  longwave  radiation,  TMP_RV  is  the 
visible  irradiance,  M  is  the  total  number  of  grid  boxes,  and  Mh  is  the  number  of  grid  box  within  the 
canopy. 

Equation  (CT.4)  is  coded  in  the  subroutine 

SV_TMP_SLP(delta,es, thetav) 

where  delta  stands  for  the  slope  ( S)  of  saturation  vapor  pressure  vs.  temperature,  es  is  the  saturation  vapor 
pressure,  and  thetav  is  the  air  temperature  at  a  specified  grid  box. 

The  vapor  pressure  (equation  (CT.5))  is  coded  in  the  subroutine 

VAPOR_PRE(e,q) 

where  e  is  the  output  vapor  pressure  and  q  is  the  specific  humidity. 

The  saturation  vapor  pressure  (equation  (CT.6))  is  coded  in  the  subroutine 
SAT_VAPOR_PRE(es,  thetav) 


Reference 

Jones,  H.  G.:  1992,  Plants  and  Microclimate  -  A  Quantitative  Approach  to  Environmental  Plant 
Physiology ,  Cambridge  University  Press,  New  York,  428  pp. 
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Norman,  J.  M.:  1979,  ‘Modeling  the  complete  crop  canopy’,  In:  Modification  of  the  Aerial  Environment  of 
Crops ,  ASAE  Monograph ,  249-3 1 1 . 
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The  Governing  Equations  for  the  Air  (GE) 


Algorithm 


Under  the  assumption  of  horizontal  homogeneity,  the  governing  equations  of  momentum,  heat,  and 
moisture  of  the  air,  after  volumetric  and  temporal  averaging,  are  given  as  (Wilson  and  Shaw,  1977; 
Raupach  and  Shaw,  1982;  Raupach  et  al.,  1985;  Gro(3,  1993;  Kaimal  and  Finnigan,  1994) 


<*»>  Aj-s  i-W  U  ,.m  4(m’v» 

(GE.l) 

(GE.2) 

dda)  LH  +  SH  4iw,d'}) 

(GE.3) 

&  PaCp,ate  & 

> 

1 

II 

"lS^ 

(GE.4) 

dt  Pa^tiZ  dz 

where  w  and  v  are  the  two  horizontal  components  of  winds,  0a  and  q  are  the  potential  temperature  and  the 
specific  humidity,  respectively,  of  the  air,  Ug  and  Vg  are  the  two  geostrophic  wind  components  representing 

the  mean  horizontal  pressure  gradients,  M  is  the  mean  wind  speed  ( =  ),  Cj)  is  the  drag 


coefficient  of  the  canopy,  and  fc  is  the  Coriolis  parameter  (  =  104  S  1 ).  The  definitions  of  the  rest  of  the 
variables  can  be  referred  to  the  section  “ Canopy  Temperature ”  In  the  above  equations,  the  brackets  and 
the  overbars  indicate  the  volumetric  and  the  temporal  means,  respectively. 

In  momentum  equations  (GE.1,2),  the  first  terms  on  the  right-hand  sides  represent  the  deviation  of  the 
flows  from  the  mean  geostrophic  winds.  It  may  be  neglected  within  the  roughness  sublayer,  since  both 
Coriolis  effect  and  mean  pressure  gradient  are  several  orders  smaller  than  the  rest  of  the  terms  for  the 
microscale  systems.  The  mean  geostrophic  winds  were  defined  using  the  mean  horizontal  pressure 
gradients  (Stull,  1994): 


u 


8 


1  dP 

fcPa  % 


1  8P 

and  v  =  +  — r-= — 7T 

*  fcPa  & 


(GE.5) 
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Notice  that  LH  and  SH  in  equation  (GE.3)  are  just  the  counterparts  in  the  governing  equation  for  the 
canopy  temperature.  This  implies  that  the  loss/gain  of  energy  on  the  foliage  due  to  latent  or  sensible  heat 


serves  as  the  source/sink  to  the  surrounding  air  within  the  canopy. 

The  last  term  of  equations  (GET -4)  represents  the  net  transport  of  the  quantity  by  the  turbulence.  This 
term  is  parameterized  using  Stull’s  transilient  turbulence  theory  ,  as  discussed  in  the  section 
“ Parameterization  of  Turbulent  Transport 

By  leaving  out  the  turbulent  transport  terms  for  separate  consideration,  the  discretization  of  equations 
(GET -4)  is  written  in  a  generic  form  as 

{Z*(t  +  At))  =  {i(t))  +  F(At,  (GE.6) 


where  £  represents  any  passive  tracer  or  quantity  (w,  v,  0,  or  q),  represents  the  external  forcings  as 


+  /[<v>-(v,)] 

-  j  CDaM(u  ) 

for  £  =  u 

CDaM(v ) 

for  %  =  v 

‘  LH+SH 

for  £  =  6 

for  £  =  q 

pucPfl^ 

LH 

PA Az 

The  asterisk  in  equation  (GE.6)  represents  the  intermediate  step  in  the  numerical  integration. 


(GE.7) 


Program 

The  subroutines  for  the  intermediate  integration  of  momentum,  heat,  and  moisture  equation  (GE.6), 
respectively,  are  named  as 

FORCE_VEL(u,v,ug,vg,M) 

FORCE_TMP(u,  v,  thetav,  thetac,  q,  SH,  LH,  RN,  tmp  a,  flxRs, 
abs__Rs,  TMP  RL,  TMP_RV,  M,  Mh) 

FORCE_HMD(q,LH,M,Mh) 

where  u  and  v  represent  the  two  components  of  winds,  ug  and  vg  represent  the  two  components  of 
geostrophic  winds,  M  represents  the  total  number  of  grid  boxes,  thetav  and  thetac  represent  the 
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temperatures  of  air  and  leaves,  respectively,  q  represents  the  specific  heat  of  air,  LH  and  SH  represent  the 
latent  and  sensible  heat,  respectively,  RN  and  RV  represent  the  net  and  visible  irradiance  within  each  grid 
box,  and  Mh  represents  the  number  of  grid  boxes  within  the  canopy. 


Reference 

GroP,  G.:  1 993,  Numerical  Simulation  of  Canopy  Flows,  Springer- Verlag,  New  Y  ork,  1 67  pp. 

Kaimal,  J.  C.  and  Finnigan,  J.  L.:  1994,  Atmospheric  Boundary  Layer  Flows  -  Their  Structure  and 
Measurement,  Oxford  University  Press,  New  York,  289  pp. 

Raupach,  M.  R.,  Coppin,  P.  A.,  and  Legg,  B.  J.:  1986,  ‘Experiments  on  scalar  dispersion  within  a  model 
plant  canopy.  Part  I.  The  turbulence  structure’,  Bound.-Layer  Meteor.  35, 21-52. 

Raupach,  M.  R.,  and  Shaw,  R.  H.:  1982,  ‘Averaging  procedures  for  flow  within  vegetation  canopies’, 
Bound.-Layer  Meteor.  22,  79-90. 

Wilson,  N.  R.,  and  Shaw,  R.  H.:  1977,  ‘A  higher  order  closure  model  for  canopy  flow’,  J.  Appl.  Meteorol. 
42,  1197-1205. 
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Parameterization  of  Turbulent  Transport  (TT) 

Algorithm 

The  translient  turbulence  theory  (TTT)  (Stull,  1984,  1985,  1986,  1990,  1993;  Stull  and  Hasegawa,  1984; 
Stull  and  Driedonks,  1987)  was  applied  to  the  parameterization  of  the  turbulent  transport  terms  in 
equations  (GE.1-4). 

Using  index  /'  and  j  represent  the  index  of  grid  box,  TTT  suggests  that  the  turbulence  mechanisms 
respond  to  the  instability  voluntarily  by  mixing  the  corresponding  quantities  between  grid  boxy  and  any 
other  grid  box  i  during  the  same  time  interval  At.  If  cy  represents  the  fraction  of  air  in  box  i  that  came  from 
j  during  time  interval  At,  then  the  new  concentration  at  box  i  after  the  mixing  can  be  calculated  by  the 
summation  of  the  mixing  from  over  all  n  grid  boxes  in  the  column,  which  is 

«('  +  A')=i>9M9£'('  +  A/)  (TT.l) 

j= 1 

which  tells  us  that  when  air  is  mixed  into  box  i  from  box  j,  the  air  carries  with  it  an  amount  of  Cy  {t.  A/) 
of  the  tracer  with  concentration  £ . ,  which  can  be  recognized  as  an  n  x  n  matrix  of  mixing  coefficient 
and  is  called  a  transilient  matrix..  The  diagonal  elements  of  the  matrix,  Cu(t,  At} ,  represents  the  fraction 

of  air  that  remains  in  box  i.  Since  ^  may  be  seen  as  an  n  x  1  matrix,  equation  (TT.l)  prescribes  simple 
matrix  multiplication. 

As  suggested  by  Stull  and  Driedonks  (1987)  followed  by  other  researchers  (Chrobok  et  al .,  1992; 
Cuxart  and  Soler,  1994),  the  following  parameterization  for  the  exchange  of  matrix  Cy  was  incorporated 
into  the  model 

mjYjWYi  for  i*j 

cv=<\-YJciJ  for/ =  y  (TT2) 

j*i 

with  the  constraints: 
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(TT.3) 


Vi  V  m>  1 

ZA=1>  2.-^  =  ! 

>1  /=1  mj 

for  the  conservation  of  air  mass  and  state,  where  Yg  is  called  the  mixing  potential,  and  zw/  and  mj  represent 
the  air  mass  in  rth  and  yth  grid  boxes,  respectively.  With  the  presence  of  the  canopy,  the  mixing  potential 
may  be  modified  (Stull,  1993)  as 

i- 1 


ij-iyno-r.) 


(TT.4) 


A=/ 


for  j  >  i ,  If  j  <1 ,  then  the  product  will  be  from  k  =  j  to  A:  =  i  - 1 ,  where  yfc  is  defined  as  the 
interference  coefficient  that  varies  between  1.0  for  total  interference  by  plants,  to  0.0  for  no  interference  of 
vertical  eddy  motions,  and  Y~  has  the  same  definition  as  Yg  (see  equation  (TT.8)  below).  The  interference 
coefficient  may  be  associated  with  the  cumulative  leaf  area  index  as 


Yu  =' 


Lk+ 1/2  A-l/2 


LAI 


(TT.5) 


The  evaluation  of  Yg  was  done  following  Stull  and  Driedonks  (1987)  using  nonlocal  form  of  TKE 
equation,  except  that  the  TKE  equation  was  defined  in  the  roughness  sublayer.  According  to  Raupach  et  al. 
(1986),  the  TKE  equation  is  given  as 

0 

j\ 

dt  dz  dz  9  dz  2 

The  nonlocal  analogous  expression  may  be  written  as 


cE  ——^7  — - -dv  g—r~,  ^{W  E)  1  „  / _ 3  _3\ 

—  =  -u'w - vw —  +  4-w - - - - —  Cna  [u  +v\-  £ 

zt  zb  zb  a  Zb  o  u  \  / 


(TT.6) 
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(TT.7) 


A  ,t 


where  0 ijy  and  V~  are  the  mean  potential  temperature,  leaf  area  density,  latitudinal  and 

longitudinal  velocity  components  between  layers  /  and  /,  respectively. 
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The  mixing  potential  is  then  prescribed  as 


W  J(A u)l  +[(AM),|(AM),|  +  (Av),|(Av),|]  +  ^CD^.M;r0 


Ml 


max 


DA,t 

for  i  *  j 

T, 

„r,„)+y« 

for  i  =  j 

(TT.8) 


and  ||  denotes  the  maximum  norm  such  as 
..  (  » 


\\Y\\  =  max 


Xr,. 

Vy= i  J 


(TT.9) 


In  the  above  formulations,  T0  represents  a  time  scale  of  turbulence,  Rc  is  analogous  to  a  critical 
Richardson  number  above  which  turbulent  mixing  is  zero  (Yjj  =  0),  D  is  a  dimensionless  factor  that  scales 
the  dissipation,  Yref  is  a  reference  mixing  potential  that  accounts  for  the  internal  mixing  within  a  grid  box, 

and  My  =  +  vjj  is  the  mean  wind  speed  between  rth  and y'th  box.. 


Under  the  circumstance  when  i  =j,  Yu  may  be  interpreted  as  the  subgrid  (internal)  mixing  potential 
for  eddies  smaller  than  the  size  of  one  grid  box.  As  explained  by  Stull  (1984)  and  Stull  and  Driedonks 
(1987),  Yu  ought  to  be  made  larger  than  any  other  element  in  the  (y)  matrix  and  the  values  of  Yy  elements 
increase  monotonically  from  the  value  of  the  upper  right-most  element  in  the  (y)  matrix  toward  the  values 
on  the  main  diagonal,  i.e.,  Yfj  >  Yf  i+l  in  order  to  avoid  the  situation  of  convective  overturn.  As  the  first 


trial,  values  of  R0  =  0.21,  D  =  1,  Tref  =  1000,  T0  =  1000  from  the  literature  were  used,  although  the 
parameters  in  the  transilient  matrix  have  to  be  determined  accordingly  by  the  sensitivity  study. 

In  equation  (TT.4),  two  terms  associated  with  the  presence  of  the  canopy  are  turbulent  transport, 


4 WE )  i 

7  ,  and  the  wake  production,  —  CDa\U  +  V  J  .  The  turbulent  transport  term  acts  as  a  significant 


sink  for  turbulent  kinetic  energy  at  the  top  of  the  canopy  and  up  to  about  z  =  2H ,  and  source  in  the 
canopy  (Kaimal  and  Finnigan,  1994;  Leclerc  et  al. ,  1990;  Meyers  and  Baldocchi,  1991),  where  both  shear 
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and  shear  production  have  fallen  to  insignificantly  low  values.  This  term  represents  the  non-local 
influences  of  TKE  and  is  important  in  maintaining  turbulent  kinetic  energy  levels  in  the  canopy.  It  also 
associates  with  the  large  values  of  skewness  observed  in  the  lower  canopy.  The  wake  production  term  is 
the  representation  of  the  drag  by  the  canopy  which  converts  mean  kinetic  energy  to  turbulent  energy.  It 
accounts  for  the  local  variations  in  shear  stress  doing  work  against  local  variations  in  mean  strain  rates.  As 
has  been  discussed  by  Kaimal  and  Finnigan  (1994),  this  term  reaches  its  maximum  in  the  upper  canopy  and 
disappears  above  the  canopy  and  toward  the  ground  surface. 

The  fluxes  of  variable  £is  given  following  Stull  (1994)  as 

V?{k)=Vt'(k-  1)  +  (^)Lc*/(4 -l)  (TT.10) 


Program 

The  subroutine  to  calculate  the  mixing  potential  (equations  (TT.6,7))  are  named  as 

MIXING_POT(Y,u,v,thetav,M) 

where  Y  represents  the  mixing  potential  as  the  output,  the  rest  of  the  variables  are  the  same  as  in  the 
section  of  u  The  Governing  Equations  for  the  Air”. 

The  subroutine  for  the  transilient  matrix  is  named  as 
TranMTRX(c,Y,M) 

where  c  is  the  transilient  matrix  as  the  output,  Y  and  M  are  the  inputs. 

The  mixing  processes  (equation  (TT.l))  were  achieved  using  subroutine 

TT_MIXING(^temp,c,M) 

where  £  represents  any  variable  (u,  v,  q,  or  q),  and  temp  is  a  temporary  variable  which  has  the  same 
dimension  as  £ 

The  fluxes  (equation  (TT.8))  were  calculated  using  subroutine 
FLUXES(£w£c,M) 

where  w£  represents  the  output  vertical  fluxes  of  the  variable  £ 
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